Inference interractions
Note
Please, click on this link to get back to the Home Page.
Here I will compare results of packages HMSC, gllvm, ecoCopula and EMtree.
Packages & data
# Packages
## Models
library(Hmsc)
library(ecoCopula)
library(PLNmodels)
library(EMtree)
library(gllvm)
## Graph
library(ggraph)
library(tidygraph)
## Plots
library(lattice)
library(corrplot)
library(gclus)
library(viridis)
# Data
data("spider")
sp_names <- colnames(spider$abund)
abund <- as.matrix(spider$abund)
env_sc <- data.frame(scale(spider$x))
## seed
set.seed(42)
## Number of chains for MCMC
nChains <- 3
## Number of cores for parallel computing
if (parallel::detectCores() > 2 & parallel::detectCores() > nChains) {
nParallel <- nChains
} else {
nParallel <- 1
}HMSC
First, let’s try HMSC. We don’t have informations about sampling design, let’s assume that each site have been visited once.
studyDesign <- data.frame(sample = as.factor(1:nrow(abund)))
rL <- HmscRandomLevel(units = studyDesign$sample)Now, we will assume that there is no interactions between covariables.
model_hmsc <- Hmsc(Y = abund, XFormula = ~soil.dry + bare.sand + fallen.leaves +
moss + herb.layer + reflection, XData = env_sc, studyDesign = studyDesign, ranLevels = list(sample = rL))Then we set up the MCMC chains.
if (.Platform$OS.type != "windows") {
thin <- 10
samples <- 500
transient <- 5 * thin
verbose <- 1
} else {
thin <- 10
samples <- 5500
transient <- 50 * thin
verbose <- 1
}And run the model.
model_hmsc <- sampleMcmc(model_hmsc, thin = thin, samples = samples, transient = transient,
nChains = nChains, verbose = 1, nParallel = nParallel)Then we diagnostic the model.
model_hmsc_post <- convertToCodaObject(model_hmsc)
par(mfrow = c(2, 2))
hist(effectiveSize(model_hmsc_post$Beta), main = "ess(beta)")
hist(gelman.diag(model_hmsc_post$Beta, multivariate = FALSE)$psrf, main = "psrf(beta)")
hist(effectiveSize(model_hmsc_post$Omega[[1]]), main = "ess(omega)")
hist(gelman.diag(model_hmsc_post$Omega[[1]], multivariate = FALSE)$psrf, main = "psrf(omega)")Not so great for chain convergence…
Let’s see the correlation with the environnemtal variables.
postBeta = getPostEstimate(model_hmsc, parName = "Beta")
plotBeta(model_hmsc, post = postBeta, param = "Support", supportLevel = 0.95)There is no correlation with the species and the bare sediment and reflexion.
OmegaCor <- computeAssociations(model_hmsc)
supportLevel <- 0.95
toPlot <- ((OmegaCor[[1]]$support > supportLevel) + (OmegaCor[[1]]$support < (1 -
supportLevel)) > 0) * OmegaCor[[1]]$mean
corrplot(toPlot, method = "color", col = colorRampPalette(c("blue", "white", "red"))(200),
title = paste("random effect level:", model_hmsc$rLNames[1]), mar = c(0, 0, 1,
0))So there is no inter-specific interactions.
gllvm
First, let’s choose the number of latent variables.
criterias <- NULL
for (i in 0:5) {
fiti <- gllvm(abund, env_sc, family = "negative.binomial", num.lv = i, formula = abund ~
soil.dry + bare.sand + fallen.leaves + moss + herb.layer + reflection)
criterias[i + 1] <- summary(fiti)$AICc
names(criterias)[i + 1] = i
}
print(criterias) 0 1 2 3 4 5
1137.907 1131.179 1170.915 1200.526 1229.728 1258.644
Based on \(AIC\), We will take 1 latent variable. Now we are fitting the model.
model_gllvm <- gllvm(abund, env_sc, family = "negative.binomial", num.lv = 1, formula = abund ~
soil.dry + bare.sand + fallen.leaves + moss + herb.layer + reflection)Let’s see some diagnostic plots.
Effect of environnemtal variables.
Species interactions
cr0 <- getResidualCor(model_gllvm)
corrplot(cr0[order.single(cr0), order.single(cr0)], diag = FALSE, type = "lower",
method = "square", tl.cex = 0.8, tl.srt = 45, tl.col = "red")Besides one specie (Arctperi), all have negative residual correlation.
ecoCopula
GCGM_mod <- manyglm(mvabund(abund) ~ as.matrix(env_sc))
GCGM_graph <- cgr(GCGM_mod)
plot(GCGM_graph, pad = 1)It’s making a small graph, some species do not seems to interact.
EMTree
We’ll start with the same model
Initialization...
Adjusting a PLN model with full covariance model
Post-treatments...
DONE!
Likelihoods: 17.61764 , 17.6573 , 17.65734 ,
Convergence took 0.27 secs and 3 iterations.
Likelihood difference = 4.241987e-05
Betas difference = 8.172318e-08
List of 5
$ edges_prob : num [1:12, 1:12] 0 0.2777 0.1111 0.1336 0.0908 ...
$ edges_weight: num [1:12, 1:12] 0 0.00761 0.00756 0.00757 0.00755 ...
$ logpY : num [1:3] 17.6 17.7 17.7
$ maxIter : num 3
$ timeEM : 'difftime' num 0.26788592338562
..- attr(*, "units")= chr "secs"
edges_prob <- EMT_model_res$edges_prob
edges_prob[edges_prob < 2/p] <- 0 # Remove links too improbable
EMtree::draw_network(edges_prob, nodes_label = sp_names, pal = "dodgerblue3", layout = "nicely",
curv = 0.1)$G
$graph_data
# A tbl_graph: 12 nodes and 20 edges
#
# An undirected simple graph with 2 components
#
# Node Data: 12 x 8 (active)
btw bool_btw bool_deg deg title name label finalcolor
<dbl> <lgl> <lgl> <dbl> <chr> <chr> <chr> <lgl>
1 8 FALSE TRUE 4 "" Alopacce "Alopacce" FALSE
2 0 FALSE TRUE 4 "" Alopcune "Alopcune" FALSE
3 0 FALSE TRUE 2 "" Alopfabr "Alopfabr" FALSE
4 0 FALSE TRUE 2 "" Arctlute "Arctlute" FALSE
5 0 FALSE FALSE 0 "" Arctperi "" FALSE
6 0 FALSE TRUE 2 "" Auloalbi "Auloalbi" FALSE
# … with 6 more rows
#
# Edge Data: 20 x 6
from to weight btw.weights neibs title
<int> <int> <dbl> <dbl> <lgl> <chr>
1 1 2 0.278 1.53 FALSE ""
2 1 8 0.369 1.31 FALSE ""
3 1 9 0.362 1.32 TRUE ""
# … with 17 more rows